Density dependence of the electronic supershells in the homogeneous jellium model 
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We present the results of self-consistent calculations of the electronic shell and supershell structure 
for clusters having up to 6000 valence electrons. The ionic background is described in terms of a 
homogeneous jellium. The calculations were performed for a series of different electron densities, 
resembling Cs, Rb, K, Na, Li, Au, Cu, Tl, In, Ga, and Al, respectively. By analyzing the occupation 
of the energy levels at the Fermi energy as a function of cluster size, we show how the shell and 
supershell structure for a given density arises from the specific arrangement of energy levels. We 
investigate the electronic shells and supershells obtained for different electron densities. Using a 
scaling argument, we find a surprisingly simple dependence of the position of the supernodes on the 
electron density. 
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I. INTRODUCTION 
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Following the observation of electronic shells.] and the 
prediction of supershells ,□ the electronic steucture of 
metal clusters has attracted much interest J-rt- Treating 
the valence electrons as an ideal Fermi gas moving in 
a suitably chosen model potential, the observed struc- 
tures can be well understoodBoEl Such an effective one- 
electron model even lends itself to a semiclassical treat- 
ment. Using a periodic orbit expansionQ_j for the density 
of states, the supershell structure can be interpreted as 
the beating pattern of the contribution of two periodic 
orbits.oH-l But obviously some arbitrariness is involved 
in choosing a "suitable" model potential to describe a 
cluster. A systematic framework for determining such 
effective one-electapn-potentials is provided by density- 
functional theory.l_!_T__l Given the arrangement of the ions 
(Born-Oppcnhcimer approximation), the effective one- 
electron potential and the total energy of a cluster can 
be determined using the Kohn-Sham formalism. The mi- 
croscopic structure of clusters is, however, not known in 
general and the determination of the ground-state geom- 
etry from first principles is a major undertaking, which 
can only be perfacmed for clusters made of not more than 
some ten atoms .EJlJ 

Therefore drastic simplifications are called for to make 
self-consistent calculations for very large clusters feasible. 
In parallel, to the early electronic structure calculations 
for solidsE-- and surfaces,--- one can, in a first approxima- 
tion, describe the ionic background by a smooth charge 
distribution, i.e. by a jellium. In the case of clusters the 
main purpose of the jellium is to restrict the valence elec- 
trons to a finite volume, that is to define the cluster size. 
Imagining a cluster to be a small drop it seems reasonable 
to choose a spherically symmetric, homogeneous jellium. 
This leads,|to the homogeneous, spherical jellium model 

(HSJM)Ma 

It turns out that assuming sphericity of the cluster 
will essentially influence only the amplitude of the shell 
oscillations. Spherical models tend to overestimate the 



amplitude of E since there are no degrees of freedom al- 
lowing fat static Jahn- Teller deformations in open-shell 
systems. E3 Such effects give rise to a fine structure be- 
tween the major shell minima. The latter, however, be- 
ing related to spherical closed-shell clusters, are properly 
described by the model. Furthermore, it was recently 
shown that even clusters with a rough surface cap_.be 
well described by a spherically averaged potential.__J It 
thus appears that assuming spherical symmetry of the 
clusters for describing the electronic shell and supershell 
structure is justified. 

The homogeneity of the jellium, however, may well be 
questioned. It has turned out that the shape of the jel- 
lium edge or, more generally, the shape of the Kohn- 
Sham potential at the surface seems to influence the su- 
pershells strongly, .sszhile leaving the electronic shells rel- 
atively unaffected. E3 In particular for high-density met- 
als large discrepancies between the predictions of the 
HSJM and the observed supershell structure have been 
observed, most notably in the case of gallium.l__l Clearly 
the simple HSJM has to be refined, e.g. by choosing a 
suitable density profi le at the surface and by including 
pseudopotentials .E_r__- In spite of that, the homogeneous 
jellium model is still used as point of reference to com- 
pare to the results of refined models. More importantly, 
a good understanding of this basic model is needed to be 
able to find the relevant features that should be included 
in an improved model, without introducing unnecessary 
complications. 

The reason for including the noble metals Au and Cu 
is to bridge the gap between the low density alkali metals 
and the high-density trivalent materials. It is clear that 
the presence of d states will affect the electronic shells. In 
the present work the results for these metals mainly serve 
for uncovering the trends in the electronic supershells as 
the electron density is varied. 

The purpose of the present paper is threefold, which is 
reflected in its organization. In Sec. ||, we briefly review 
the homogeneous jellium model and present the results of 
our calculations. These include the oscillating part E(N) 
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of the total energy for densities resembling Cs, Rb, K, Na, 
Li, Au, Cu, Tl, In, Ga, and Al, respectively, for clusters 
having 125 up to 6000 valence electrons. The position of 
the minima in E (magic numbers) are listed explicitly. 
Furthermore, a method for determining the location of 
the nodes in the supershell oscillation is given and the 
corresponding results for the clusters mentioned above 
are shown. This collection of data could serve as refer- 
ence for jellium results— .pjtandipgj previously published 
results considerably.OEaila'LJ~El'E3 

Focusing on the occupation of energy levels as the clus- 
ters increase in size, a simple interpretation of the origin 
of supershell structure can be given. This is done in 
Sec. IIL Analyzing the states at the Fermi energy, we 



show that the supershells arise from the fact that the 
difference in angular momentum A^ for orbitals of given 
number n of nodes in the radial wave function increases 
with cluster size. Since only integer values of Al corre- 
spond to physical degeneracies, there is a periodic change 
from degenerate to nondegenerate energy levels as the 
clusters grow larger. We give a semiclassical explana- 
tion of this observation, which provides an alternative to 
arguments using periodic orbit expansions. 

In Sec. IV we analyze how the shell and supershell 
structure depends on the electron density. We show that 
for a homogeneous jellium edge the softness of the po- 
tential at the surface is fairly independent of the density. 
It thus turns out that only the Wigner-Seitz radius r s 
(related to the electron density by 3/(47rrg)) is the rele- 
vant length scale for the electronic structure. We actually 
find that the supershells are linearly shifted with l/r s , 
while the magic numbers are practically independent of 
the electron density. This suggests a general mechanism 
of how the surface softness influences the supershell struc- 
ture. Finally, a fit to our data allows the results of jellium 
calculations to be estimated without actually having to 
perform the lengthy computations. 

Throughout the paper we give lengths in Bohr radii 
(ao) and energies in Rydbergs (Ry)- 



To find the total energy E(N) of a jellium cluster hav- 
ing N valence electrons we solve the Kohn-Sham equation 



( - A + V KS (r ) = e^(f ) 



(1) 



sclf-consistently. The Kohn-Sham potential is given by 
the sum of the electrostatic potentials arising from the 
jellium and the electronic charge density, respectively, 
and the exchange-correlation potential 

V KS {r) = V z (r) + 2 / dV ^Kjt + V xc (f). 

J \r-r'\ 

The electron density n e i(f) is determined from the N 
lowest-lying eigenfunctions of (fy) 



n e i(r ) = > W„(r 

fi:occ 



Having reached self-consistency, the total energy is given 
by 



E(N) =Y. e »- J d 3 rV K s{r)n e i(f) 



+E xc [n] + Ecoui[n], 



(2) 



where E xc is the exchange-correlation energy functional 
and Ecoui is the electrostatic self-energy of the total 
(both jellium and electronic) charge density. 

To treat large clusters more efficiently, instead of ac- 
tually finding all the eigenfunctions of the Kohn-Sham 
equation, we determine the electron density n e i(r) and 
the sum ^2 of the one-electron energies by contour in- 
tegration from the one-electron Green's function of (Q). 
Exploiting spherical symmetry, the observables are given 
by 



n e i(r) = 



E 2 -&£if9l(E;r,r)dE (3) 

Ep 



II. JELLIUM CALCULATIONS 
A. Method 

In the homogeneous jellium model a given material 
is characterized by the average density ribuik of the va- 
lence electrons. Later on, it will be convenient to use 
the Wigner-Seitz radius r s , in terms of which the bulk 
electron density is given by n^ik = 3/(47rr^). To make 
the description of the clusters independent of the valence 
of the metal under consideration, we will use the num- 
ber N of valence electrons, not the number of atoms, to 
denote the size of a cluster. In the spherical, homoge- 
neous jellium model the jellium density is then given by 
nj(r) — ribuik O(i?o— r), where Rq = r s A 1 / 3 is the radius 
of the cluster. 



oo 

E^-^/drr 2 ^ I gi ( E ;r,r)EdE. (4) 

" ' E F 

Here gi(E;r,r) is the trace of the Green's function for 
the radial Kohn-Sham equation with angular momentum 
I. It is constructed from the two solutions uf(E;r) and 
uf{E\r) which are regular at the center and at infinity, 
respectively 



u<(E;r) u>(E;r) 
W{E) 



where W{E) = uf{E;r) u>{E]r)~u<{E;r) uf'{E;r) is 
the Wronskian of (|l|). 

At zero temperature the integration contour would 
have to go through the Fermi energy Ep and we would 
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have to deal with poles on the real axis. We therefore in- 
troduce a fictitious finite but small temperature and mul- 
tiply the Green's function in Eqs. (|J) and (||) by a Fermi- 
Dirac distribution. The contour can then be closed well 
above Ep, where the poles on the real axis have negligi- 
ble weight. The Fermi-Dirac distribution, however, has 
additional poles at Ep + (2n + l)nikT, for integer n. 
Their contributions have to be subtracted from the con- 
tour integral. Working with fractional occupations has 
the further advantage of reducing oscillations in the oc- 
cupation numbers during the iterative process of reaching 
self-consistency. 

The main reason for using the Green's function ap- 
proach to solving the Kohn-Sham equation is its com- 
putational efficiency for large clusters. Since the radial 
mesh grows with the cluster radius i?o = r s TV 1 / 3 , the 
computation of the Green's function scales with A 1 / 3 . To 
determine n e i and we have to perform l m ax oc Ro 

integrations, for which the integration contour is inde- 
pendent of N. Hence the computer time grows as TV 2 / 3 . 
This has to be compared to the traditional approach, 
which for each I involves finding n max oc R eigenstates, 
thus resulting in a computational complexity of O(N). 

Because of the the large density gradients at the clus- 
ter surface, it seems highly desirable to go beyond local- 
density approximation to the exchange-correlation en- 
ergy. Fig. [j] shows the self-consistent Kohn-Sham po- 
tentials for a Na5ooo jellium obtained using local density 
approximation and two-different generalized gradient ap- 
proximations (GGA's).E3 : EJ Unfortunately, it turns out 
that the potentials found using a GGA exhibit an un- 
physical shoulder at the surface- We therefore stick to 
the local-density approximation^ for all subsequent cal- 
culations. 



B. Results 

We have performed self-consistent calculations for jel- 
lium clusters with electron densities corresponding 
to Cs, Rb, K, Na, Li, Au, Cu, Tl, In, Ga, and Al (see 
Table |). Cluster sizes range from 125 to 6000 valence 
electrons for each density considered. The calculations 
yield the total energy E{N) as a function of cluster size. 

Neglecting the discreteness of energy levels (i.e., disre- 
garding shell effects) the total energy is approximated by 
the asymptotic expansion 



E{N) =a 1 N + a 2 N 2/3 + a 3 N 1/3 + 



(5) 



In this Thomas-Fermi-like expansion the leading coeffi- 
cient a\ resembles the bulk energy of the system— The 
next coefficient a 2 is related to the surface energy.El The 
electronic shells and supershells reflect the deviations of 
the total energy E(N) from this smooth function. These 
are given by the oscillating part of the total energy, de- 
fined by 
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FIG. 1. Self-consistent Kohn-Sham potentials for a ho- 
mogeneous, spherical Nasooo jellium (Vks = V CO ui + V xc )- 
The results shown here differ only in the way the ex- 
change-correlation term was treated. The potential displayed 
in the first plot (LDA) was obtained using the local-density 
approximation in jfeke parametrization given by Vosko, Wilk, 
and Nusair (Ref.B). The second (GGA-PW) and third 
(GGA-EV) plot show the potentials obtained using gener- 
alized gn&dient approximations as givep-iby Perdew-Wang et 
al. (Ref.Ej) and Engel and Vosko (RefEU), respectively. 



Thus, to extract E(N) from the total energy E{N) 
obtained from our self-consistent calculations we need to 
know the smooth part E{N) of the total energy. To this 
end we have determined the parameters dj in the asymp- 
totic expansion (||) by a least-squares fit. To check the 
quality of our fit, we compare the fit parameters to known 
properties of jellium systems. The leading parameter a\ 
should be equal to the bulk energy per electron of a ho- 
mogeneous, neutral electron gas 



E(N) = E(N) - E(N). 



(6) 
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The comparison of the parameters a\ from the fit to tbuik 
is shown in Fig. |^. 
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FIG. 2. Leading term in the asymptotic expansion of the 
total energy E(N) extracted from our calculations for ho- 
mogeneous jellium spheres of densities n e i — 3/(47rrJ). For 
comparison the dotted line shows the energy of a homoge- 
neous electron gas in local density approximation usins-the 
parametrization given by Vosko, Wilk, and Nusair (RefE3). 

The next to leading term a-i TV 2 / 3 should resemble the 
total energy of the jellium surface, hence 02 should be 
related to the surface energy a by 



Q2 

47rr? 



(7) 



Fig. y shows the comparison between the surface en- 
ergy obtained from the fit parameter 02 and the surface 
energies for jellium, surfaces calculated by Lang-KohntH 
and Perdew et alJB. The discrepancies between the sur- 
face energies found by Lang-Kohn and the other data 
are probably due to the use of a different exchange- 
correlation energy functional: Lang and Kohn used the 
Wigner approximation to the exchange-correlation en- 
ergy, while the treatment of exchange and correlation 
in the work of Perdew et al. and in our calculations was 
based on the results of Mpiite Carlo simulations of the 
homogeneous electron gas a 

Having found the parameters of the asymptotic ex- 
pansion (|5j) we can readily determine the oscillating part 
E(N) of the total energy. The results from our self- 
consistent calculations for the densities listed in Table 
U are shown in Fig. [|. Plotting E over N 1 / 3 , which is 
proportional to the cluster radius, we find regular oscil- 
lations (shell structure). The amplitude of these oscilla- 
tions is modulated, resembling a beating pattern (super- 
shell structure). 

The local minima in the oscillating part E of the to- 
tal energy correspond to clusters that are expected to be 
exceptionally stable. We call the number of valence elec- 



trons TV in these clusters magic numbers. They are listed 
in Table [jj. 
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FIG. 3. Next to leading term in the asymptotic expansion 
of the total energy E(N) extracted from our calculations for 
homogeneous jellium spheres of densities n e i = 3/(47rrf ) (di- 
amonds). For comparison the surface energies a for jellium 
surfaces as found by Lang and Kohn (Ref.E-J, Table II) and 
by Perdew et al. (Ref.EJ, Table XII) are shown. The dotted 
line gives an interpolation of the Lang-Kohn data. 

It turns out that in the regions where the amplitude of 
the shell oscillations is large, the magic numbers for dif- 
ferent jellium densities are comparable. Comparing the 
sequences of corresponding minima in E, we find that 
the magic numbers grow slightly as the Wigner-Seitz ra- 
dius decreases. In the regions where the amplitude of 
the shell oscillations is small, similar sequences of magic 
numbers exist. Here, however, they do not extend over 
the whole range of jellium densities. Instead we find that 
a sequence starting at small densities ends at some inter- 
mediate density, while a new sequence, shifted by about 
half a period, evolves and extends to larger jellium den- 
sities. 

For all jellium densities we find a pronounced modu- 
lation in the amplitude of the shell oscillations. Even a 
superficial look at Fig. [| shows that the beating pattern 
in E is shifted towards larger cluster sizes as the jellium 
density increases. To characterize the supershell struc- 
ture we use the location of the minima in the amplitude 
of the shell oscillations. The position of these supernodes 
is given by minima in the envelope of E(N). The enve- 
lope is obtained eliminating the rapid shell oscillations 
from ^(iV 1 / 3 ) using a low-pass Fourier filter. Naturally 
such a filtering scheme introduces an uncertainty in iV 1 / 3 
which is of the order of a period of the shell oscillation. 

For all jellium densities we have considered, exactly 
two supernodes fall in the range of clusters having 125 
- 6000 valence electrons. Using the procedure described 
above, we have determined the position of these supern- 
odes. The results are listed in Table [II. 
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FIG. 4. Oscillating part E of the total energy for jellium 
clusters as a function of the number of valence electrons N. 
At the minima of E the corresponding number of electrons is 
indicated (magic numbers). These numbers are also listed in 
Table [n|. The positions of the supershell nodes are indicated 
by open triangles. The numbers given in parentheses below 
these are the positions of the supernodes given as N 1 ^ 3 (cf. 
Table Q. 
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Cs 


Rb 


K 


Na 


Li 


Au 


Cu 


Tl 


In 


Ga 


Al 


r s (do) 


5.62 


5.20 


4.86 


3.93 


3.25 


3.01 


2.67 


2.48 


2.41 


2.19 


2.07 


TABLE I. Wigner-Seitz radii r s for the different metals which 


were used in the self-consistent calculations to define the 


jellium density. 
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Tl 


111 


Pa 


Al 

r\l 




137 


137 


138 


138 


138 


138 


138 


138 


138 


138 


138 




186 


186 


186 


189 


196 


196 


196 


198 


198 


198 


198 




254 


254 


254 


254 


254 


254 


254 


261 


267 


267 


267 




338 


338 


338 
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338 


338 


339 


339 


339 


339 




437 


438 


438 


439 


440 


440 


440 


441 


441 


441 


441 




539 


543 


545 


551 


556 


556 


557 


558 


561 


561 


561 




617 


618 


619 


634 


640 


638 


638 


639 


639 


639 






668 


673 


674 


676 


676 


677 


686 


693 


693 


693 


693 




751 


751 


752 


756 


759 


783 




















832 


832 


832 


832 


831 


831 


834 


834 




910 


911 


912 


912 


912 


914 


921 


924 


927 


927 
















1013 


1012 


1011 


1011 


1011 


1011 




1078 


1079 


1080 


1096 


1100 


1100 


1100 


1101 


1101 


1101 


1101 






















1215 


1218 




1281 


1282 


1282 


1284 


1284 


1286 


1314 


1314 


1314 


1314 


1314 




1503 


1503 


1503 


1507 


1515 


1516 


1518 


1521 


1521 


1524 


1554 




1747 


1751 


1755 


1760 


1760 


1760 


1760 


1767 


1770 


1776 


1779 




2018 


2019 


2019 


2019 


2039 


2047 


2048 


2049 


2049 


2049 


2049 




2322 


2324 


2326 


2330 


2333 


2334 


2334 


2337 


2346 


2367 


2367 




2654 


2653 


2653 


2656 


2665 


2670 


2672 


2673 


2673 


2679 


2682 




2909 


2906 


2903 


2913 


















3070 


3063 


3056 


3025 


3028 


3028 


3028 


3030 


3030 


3042 


3048 




3261 


3267 


3271 


3279 


3283 


3288 






















3409 


3418 


3436 


3438 


3438 


3438 


3438 










3510 


















3607 


3619 


3645 


3676 


3692 


3695 




3708 


3708 


















3865 


3847 


3849 


3849 


3861 


3882 




4063 


4067 


4068 


4068 


4088 


4112 


4152 


4152 


4155 


4155 


4158 


















4332 


4320 


4317 


4323 




4544 


4554 


4561 


4569 


4570 


4570 


4570 


4569 


4569 


4593 


4656 




5015 


5025 


5035 


5068 


5095 


5097 


5101 


5106 


5106 


5109 


5109 
















5471 






5511 


5502 




5567 


5565 


5565 


5582 


5601 


5607 


5630 


5655 


5658 


5670 


5670 


TABLE II. List of the mag 


ic numbers for jellium clusters of different densities. Listed are the numbers of valence electrons 


N in the range from 125 to 6000 for which the oscillating part of the total 


energy has 


a minimum. The 


magic numbers for 


the different jellium clusters are arranged in sequences of comparable size. 


In the regions where the amplitude of the shell 


oscillations is large 


(cf. Fig. U 


) these sequences extend 


over the whole range 


of densities. In the vicinity of the nodes 


in the 


supershell structure this is no 


longer true. 


Here sequences starting at small jellium densities 


end at intermediate densities, 


while a new sequence, shifted by about half 


a period, evolves and extends to larger densities. 










Cs 


Rb 


K 


Na 


Li 


Au 


Cu 


Tl 


In 


Ga 


Al 




8.39 


8.43 


8.47 


8.95 


9.15 


9.30 


9.68 


9.79 


9.85 


10.13 


10.29 


N? 3 


14.59 


14.67 


14.75 


15.11 


15.45 


15.62 


15.86 


16.01 


16.11 


16.37 


16.49 



TABLE III. Position of the nodes in the supershell structure for jellium clusters of different densities. The position of the 
supernodes was determined by finding the minima in the envelope of E(N). N denotes the number of valence electrons in the 
cluster. 
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III. MICROSCOPIC ANALYSIS 

A cluster with N valence electrons will be more stable 
than neighboring clusters if its highest occupied orbital 
is completely filled. The increase in stability can be char- 
acterized by the gap between the highest occupied and 
the lowest unoccupied level. In general, the width of 
this gap rapidly decreases with growing cluster size. A 
set of degenerate levels will however open a larger gap 
in the spectrum. The electronic shell structure can thus 
be traced back to degeneracies in the energy spectrum 
of the cluster potential at the Fermi energy Ep. In the 
following we therefore focus on the Kohn-Sham energies 
e Ut i next to Ep. In doing so we identify the degeneracies 
that cause the strong shell oscillations and find a mech- 
anism that explains why there are no such degeneracies 
near the supernodes. 

The potentials in our self-consistent calculations are 
spherically symmetric. Therefore we can characterize 
the energy levels by two quantum numbers: the number 
n of nodes in the radial wave function and the angular 
momentum quantum number I. The degeneracy of an 
eigenenergy e n j is given by 2(21 + 1). Hence states with 
large I have a stronger influence on the shell structure 
than those with small angular momentum. 

Since in our self-consistent calculations the Kohn- 
Sham levels are populated according to a Fermi-Dirac 
distribution 



2(2/ + 1) 



1 + e (e n ,i-E F )/kT 



degenerate levels e raj ; near the Fermi energy Ep are filled 
simultaneously. By "filling" a level we mean that the oc- 
cupation of e n j changes with cluster size (i.e. with the 
total number N of valence electrons) : 

AN en>l (N) = N tn , (N) - N entl (N - 1) + 0. 

In practice, we do not consider all individual AiV £fl , , but 
only the sums 



For a given angular momentum I, the first level to be 
filled as the cluster size increases will be the one with 
the nodeless radial function: eo t i, followed by those with 
quantum numbers n = 1, 2, . . .. Since these levels are well 
separated, their filling will be indicated by small humps 
in ANi, which otherwise vanishes. Fig. [3] shows the ANi, 
shifted vertically by I, as functions of iV 1 / 3 . 

Since orbitals having the same energy are filled simul- 
taneously, finding several humps for the same cluster size 
means that the corresponding orbitals e Ut i are degener- 
ate. For the the leading levels (i.e. those with large angu- 
lar momentum) these degeneracies are indicated in Fig. ||. 
The patterns are similar for all jellium densities. 
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FIG. 5. Filling of the orbitals with increasing cluster size 
N in the self-consistent calculations for jellium clusters. ANi 
is the change in the occupation of the orbitals with angular 
momentum I. In the plot, the AiV; are shifted by I. The small 
humps in the horizontal lines indicate a change in the occupa- 
tion of an orbital with the corresponding angular momentum, 
i.e., the filling of that orbital. Since degenerate orbitals are 
filled simultaneously, the existence of several such humps for 
a given cluster size N shows that the corresponding orbitals 
are degenerate. For the orbitals with large angular momen- 
tum these degeneracies are indicated by lines connecting such 
humps. The first hump for a given I corresponds to the or- 
bital with nodeless radial wave function n = 0. Likewise, the 
orbital connected with the second hump has the quantum 
number n = 1. As can be seen from the plots, the humps for 
orbitals with given n are arranged into nearly straight lines 
of slightly different slope. 

The main difference is that the pattern is shifted towards 
larger cluster sizes as the Wigner-Seitz radius decreases. 
We find that the first maximum in the shell amplitude 
is caused by the degeneracy of the states eo,; ~ £i,z-3 
followed by eo,; ~ £2, (-6- The second maximum in the 
supershell structure is mainly due to the degeneracies 
« ea,l-3, supported by e ,z ~ ei,2-4 and ei,j « e 3: ;_ 6 . 
To find out why there are nodal regions in E, we first 
observe that the humps with a given quantum number n 
are arranged into nearly straight lines (see Fig. ||). The 
slope of these occupation lines is larger for smaller n, i.e., 
the distance Al between two such lines increases with 
cluster size N . But since the angular momentum is 
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FIG. 6. Position of the energy levels e n> i = h 2 kn i/2m in 
an infinite spherical potential well (spherical cavity) of radius 
Ro- For comparison with Fig. [|, we give kRo on the abscissa 
and the angular momentum I on the ordinate. The position 
of the energy levels is indicated by crosses. For small n, these 
can be well described by the semiclassical expansion that is 
described in the text. The classical limit L max = hkpRo is 
given by the dotted line. 

quantized, degenerate energy levels exist only for integer 
values of Al. Therefore regions where energy levels are 
degenerate are separated by regions without degenera- 
cies. Hence the supershell structure can be traced back 
to the differences in the slopes of the occupation lines. 

The origin of these lines can be understood from a sim- 
ple classical argument. Assuming that the Fermi energy 
Ep is fairly independent of the cluster size, the maximum 
angular momentum for an electron in a potential well of 
radius Rq = r s TV 1 / 3 is given by 



L m ax(N) = hk F R , 



(8) 



where Tikp — {2mEp) 1 / 2 is the Fermi momentum. Using 
semiclassical quantization, we can go beyond this crude 
approximation. The angular momentum is then written 
as L — h(l + 1/2), where I is the angular momentum 
quantum number. Introducing the parameter 
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1 + 1/2 
kpRa ' 



(9) 



which is unity in the classical case QSJ) , the Bohr- 
Sommerfeld condition for an infinitely deep potential well 
of radius Rq reads 



(l + \) (J^- 1 - arccos( 7 )J = tt (n + 3/4) . (10) 

Expanding this equation to first order around the classi- 
cal case 7 = 1, we can solve for the angular momentum 
of a state with quantum number n at the Fermi level: 



kpRf) — 



(37r(n + 3/4) 



2/3 



Fig. H shows the angular momentum quantum number 
I as a function of kpR . It turns out that, for small 
quantum numbers n, the approximation ([l|) matches the 
quantum mechanical result quite well. The deviations for 
larger n are due to the expansion around 7=1, not to 
the semiclassical approximation itself. Since a state e Wj ; 
at the Fermi level will be filled as the cluster size in- 
creases, the crosses in Fig. || can be identified with the 
humps in Fig. |5|. Hence the expansion (|ll|) provides, for 
small n, a good approximation to the occupation lines 
for the infinite potential well. 



IV. SCALING ANALYSIS 

In the present section we analyze the systematic vari- 
ations in the shell and supershell structure as the jellium 
density is changed. The key to such an analysis is the 
identification of relevant scales. Expressing all data in 
terms of these scales, i.e. dealing with dimensionless re- 
duced quantities, allows the results for different materials 
to be conveniently compared. 

Obviously, the Wigner-Seitz radius r s is such a scale. 
It is a fundamental length describing how the volume of 
the cluster depends on the jellium density. Via the rela- 
tion 




it is also related to a basic energy scale. As a simple 
example we note that the overall amplitude of the shell 
oscillations E is proportional to \jr\. This is shown in 
Fig. | 

From the shape of the self-consistent Kohn-Sham po- 
tentials we can identify a second length scale. As can be 
seen from Fig. |§|, for a given jellium density the shape 
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FIG. 7. Scaling of the amplitude of the shell oscillations 
with the electron density n e i = 3/(47rrJ). The diamonds give 
the maximum amplitude of the oscillating part E extracted 
from our self-consistent calculations. The dotted line gives 
a fit of these data points with the one-parameter function 
const. /r s . 
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FIG. 8. Self consistent Kohn-Sham potential for sodium 
and gallium jellies of different sizes. To compare the surface 
region of the potentials, we chose r minus the radius Ro of 
the corresponding jellium as the abscissa. As can be seen, the 
shape of the potential at the surface of the clusters is fairly 
independent of the size N. The horizontal lines indicate the 
Fermi energies of the corresponding systems. 

of the potentials around i?n is quite independent of clus- 
ter size. For each material it can be characterized by the 
width of the surface region, which hence is a good can- 
didate for an additional length scale. We will denote the 
surface width by a. 

We might assume then that the changes in the elec- 
tronic shells and supershells can be described in terms 
of the dimensionless parameter a/r s . For a rough esti- 
mate of how a varies with the jellium density, we have 
performed a simple Thomas- Fermi calculation for a plane 
jellium surface nj(z) = nj0(— z). Given an area A the 
surface energy is 



o[n e i{z) 



GO , 

/ (e(n e i(z))n e i(z) - e(m(z))nj(z) 



dz 



A 



where e(n) is the Thomas-Fermi energy functional in- 
cluding—the second order gradient correction to kinetic 
energy,Ej exchange energy and Hartree energy. To 
estimate a, we minimize the surface energy using a 
parametrized electron density, where the only parame- 
ter is the surface width 



n el (z) = m f(z/a). 
A good choice for the function f(z) isl 
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FIG. 9. Width a op t of the surface region as a function of 
electron density. a opt was obtained from a Thomas-Fermi 
calculations (including gradient corrections and exchange) 
for a parametrized electron density. The energy functional 
included second-order gradient corrections to the kinetic en- 
ergy, Coulomb energy, and exchange energy (full line), 
eluding also exchange in the X a approximation (Ref.t 
yields the dotted curve. 

The resulting surface width a op t as a function of the 
Wigner-Seitz radius r s is shown in Fig. ^. It is fairly 
independent of the jellium density, in agreement with 
the results of the self-consistent Kohn-Sham calculations 
shown in Fig. || Hence r s seems to be the only relevant 
length scale in the cluster problem. 

The dependence of the shell and supershell structure 
on the scaling parameter a/r s is revealed by plotting the 
magic numbers and the supernodes against l/r s . Fig. [Io| 
shows the resulting pattern, which is surprisingly simple. 
The magic numbers are practically independent of the 
jellium density, while the supernodes are linearly shifted 
with X/r s . A closer look at the supernodes reveals that 
their position can be accurately interpolated by a linear 
function. For some given jellium density a reliable es- 
timate of the location of the minima in the supershell 
structure can be read off Fig. O, without having to per- 
form a time-consuming self-consistent calculation. 

It is interesting to note that the first and the second 
supernode are shifted by the same amount. Writing the 
oscillating part of the total energy as a beating pattern 



E 



(N) = cos (/iAT 1/3 - pi) + cos (/ 2 A^ 1/3 - </? 2 ) 



2 cos 



x cos 



/ /i-/a „i/3 _ <Pi ~<P2 
V 2 2 

fl+h N l/3 _ 



(12) 



suggests a simple picture for understanding such a par- 
allel displacement. Since the supershell oscillation is de- 
scribed by the first cosine in ( fj"2| ) , the position of the k th 
supernode is given by 
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FIG. 10. Position of the shell minima (crosses) and of the 
supernodes (diamonds) as a function of l/r s . Correspond- 
ing shell minima and supernodes are connected by full and 
broken lines, respectively. 




1/r. (1/a ) 1/r, (1/a„) 



FIG. 11. Shift of the first and second supernode as a 
function of l/r a . The error bars are due to the uncertainty 
in determining the the envelope of E(N). The lines give a 
linear fit of the data, the parameters of which are displayed 
in the plot. 

^i /3 = (2k + 1)tt <£i - tf2 

fi — h h~ h 

In this picture a shift of the supernodes, which leaves the 
period of the oscillation unchanged, naturally arises from 
a variation of the phases ip\ and tp 2 . In a_ semiclassical 
approach, using a periodic orbit expansion,0 these phases 
are determined by the cluster surface. This is consistent 
with the observation that the supershell structure is sen- 
sitive to the shape of the surface potential.ErE3 Using an 
expansion of the relevant semiclassical integrals in terms 
of the scaling parameter a/r s (leptodermous expansion) 
it can be shown that the above interpretation indeed de- 
scribes the mecianism underlying the observed shift of 
the supernodes 



V. CONCLUSIONS 

We have analyzed how the oscillating part E of the to- 
tal energy extracted from self-consistent calculations de- 
pends on the electron density. It turns out that the magic 
numbers are fairly independent of the density, while the 
supershell oscillations are shifted towards larger cluster 
sizes as the Wigner-Seitz radius decreases. Focusing on 
the energy spectrum near the Fermi level Ep, we have 
found a simple pattern in which the orbitals are filled: 
the occupation lines. They can be understood using a 
semiclassical expansion. Qualitatively there is no differ- 
ence between clusters of different density. The strong 
oscillations in E are caused by degeneracies of energy 
levels close to Ep. The supershells are a consequence of 
the different slopes of the occupation lines, which imply 
transition zones between regions where levels are degen- 
erate. The shift of the supernodes with increasing density 
can be understood from an analysis of the physical scales. 
We have identified r s as the relevant length scale for the 
cluster problem. It turns out that the location of the 
supernodes is linearly shifted with l/r s . 

Furthermore, the identification of r s as the typical 
length scale for the supershells suggests_a justification of 
the ad hoc procedure proposed in Ref.B3 to improve the 
results of jellium calculations for gallium clusters. There 
it was found that the introduction of a non-homogeneous 
jellium background is essential for treating Ga^ clusters, 
while alkali-metal clusters are well described by a ho- 
mogeneous jellium. With the help of the above scaling 
argument that puzzle may be resolved. Assuming that 
the typical length scale for features in the jellium is the 
ionic radius r at , while the length scale for the electrons is 
the Wigner-Seitz radius r s , we find that the importance 
of inhomogeneities increases with the number of valence 
electrons Z va \ oc (r at /r s ) 3 . Hence, for trivalent materials 
a soft jellium surface should have a stronger influence on 
the supershells than for alkali metals. 
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